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PROBABILISTIC  MULTI-HYPOTHESIS  TRACKING 


L  INTRODUCTION 


In  a  multitarget,  multimeasurement  environment,  knowledge  of  the  measurement-to-track 
assignments  is  typically  unavailable  to  the  tracking  algorithm.  This  report  is  a  probabilistic 
approach  to  the  measurement-to-track  assignment  problem;  that  is,  measurement  assignments  are 
modeled  as  discrete  random  variables.  Measurements  are  not  assigned  to  specific  tracks  as  in 
traditional  multi-hypothesis  tracking  (MHT)  algorithms;  instead,  the  probability  that  each 
measurement  belongs  to  each  track  is  estimated  using  an  empirical  Bayesian  algorithm.  The 
probabilistic  multi-hypothesis  tracking  (PMHT)  approach  proposed  in  this  report  treats  target 
states  and  measurement  assignments  as  continuous  and  discrete  random  variables,  respectively, 
and  defines  an  appropriate  joint  density  on  these  variables.  The  PMHT  estimation  algorithm  is  in 
the  class  of  so-called  empirical  Bayesian  methods  (reference  1);  that  is,  it  is  a  hybrid  maximum 
a  posteriori  (MAP)  and  maximum  likelihood  (ML)  algorithm.  The  PMHT  estimates  are  joint 
estimates  of  target  states  and  measurement-to-track  assignment  probabilities.  This  report  expands 
earlier  work.  (See  references  2  and  3.) 

The  PMHT  algorithm  requires  neither  enumeration  of  measurement-to-track  assignments  nor 
pruning  because  all  measurements  are  assigned  to  all  tracks.  A  measurement  is  weighted  with  the 
estimated  measurement-to-track  assignment  probability  before  it  is  assigned  to  a  target  track.  The 
same  measurement  receives  different  weights  for  different  tracks.  Because  pruning  is  not 
required,  the  PMHT  algorithm  is  an  optimal  empirical  Bayesian  multitarget  tracking  algorithm 
under  idealized  assumptions. 

The  computational  complexity  of  the  PMHT  algorithm  is  much  less  than  the  exponential 
complexity  of  MHT  algorithms,  which  use  exhaustive  enumeration.  Current  MHT  algorithms 
combine  iterative  methods  for  enumeration  with  careful  scoring  and  pruning  (e.g.,  measurement 
gating  and  branch  elimination)  to  reduce  computational  complexity  and  computer  storage 
requirements  to  manageable  levels.  The  PMHT  algorithm  presented  in  this  report  is  potentially 
more  advantageous  than  the  MHT  algorithms  because  the  PMHT  algorithm  is  amenable  to 
parallel  computation  on  high-performance  computer  architectures.  Specifically,  PMHT 
measurement-to-track  assignment  probabilities  greatly  reduce  the  need  for  branching  (i.e., 
algorithm  flow  is  smooth  and  relatively  uninterrupted). 

One  of  the  key  ideas  behind  the  PMHT  algorithm  is  the  avoidance  of  “hard”  measurement-to- 
track  assignment  decisions.  It  may  be  argued  that  hard  assignments  should  be  avoided  because 
they  are  equivalent  to  statistical  decisions,  that  all  statistical  decisions  are  opportunities  for  error, 
and  that  erroneous  decisions,  even  if  rarely  committed,  necessarily  increase  estimation  error. 
However,  as  will  be  demonstrated,  the  formulation  of  the  PMHT  probability  structure  does  not 
support  this  argument  because  the  optimal  joint  MAP  estimate  of  target  states  and  measurement 
assignments  comprise  state  estimates  and  hard  assignments.  One  of  the  primary  technical 
contributions  of  the  PMHT  approach  is  that  it  avoids  the  computational  complexity  of  MAP 
estimation  (i.e.,  the  theoretical  basis  of  MHT  algorithms)  and,  yet,  is  optimal  in  the  empirical 
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Bayesian  sense.  The  change  in  the  definition  of  the  likelihood  function  that  results  from  treating 
assignments  as  discrete  random  variables  thus  has  a  dramatic  effect  on  the  resulting  estimation 
algorithm  and  its  computational  complexity. 

This  report  consists  of  ten  numbered  sections,  the  contents  of  which  are  indicated  below: 

•  Section  1  introduces  the  PMHT  algorithm. 

•  Section  2  introduces  the  PMHT  observer  and  presents  a  theoretical  overview  of  the 
PMHT  estimation  method. 

•  Section  3  states  the  likelihood  structure  of  the  PMHT  observer.  Bayesian  inference 
networks  (BINs)  are  introduced  as  a  road  map  for  representing  the  conditional  independence 
assumptions  that  are  central  to  the  methods  of  this  study.  The  validity  of  the  PMHT  approach  for 
nonlinear  non-Gaussian  systems  is  attributed  to  the  fact  that  these  conditional  independence 
assumptions,  displayed  graphically  by  the  BIN,  are  nonparametric  assumptions  and,  therefore,  are 
not  limited  to  linear-Gaussian  parameterizations. 

•  Section  4  gives  a  derivation  of  the  PMHT  algorithm  using  the  expectation-maximization 
(EM)  method.  The  derivation  assumes  familiarity  with  the  EM  method. 

•  Section  5  states  the  PMHT  algorithm  in  recursive  form.  The  linear-Gaussian  case  was  first 
presented  (without  proof)  in  reference  2. 

•  Section  6  discusses  Fisher  information  matrices  (FIMs)  and  multitarget  observability.  FIMs 
relate  to  PMHT  estimation  error,  and  the  multitarget  observability  relates  to  PMHT  algorithm 
convergence. 

•  Section  7  discusses  several  variations  and  extensions  to  the  PMHT  approach,  including 
adaptive  covariance  estimation. 

•  Section  8  compares  the  PMHT  system  with  other  multitarget  tracking  algorithms. 

•  Section  9  presents  a  PMHT  example. 

•  Section  10  provides  concluding  remarks  and  recommendations  for  further  study. 
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2.  THEORETICAL  OVERVIEW 


Section  2  uses  an  abbreviated  form  of  the  notation  used  in  the  rest  of  this  report.  A  complete 
description  of  the  issues  mentioned  in  this  section  is  presented  in  sections  3  and  4. 

Measurements  are  outputs  from  a  sensor  signal  processor,  and  the  task  of  the  postprocessor 
is  to  compute  target  tracks  from  the  available  measurements.  PMHT  is  an  algorithm  for  the 
postprocessor.  Adopting  the  observer  concept  is  useful  when  studying  the  multitarget  tracking 
problem.  The  observer  is  not  particularly  interesting  for  single-target  tracking,  but  its  value  for 
multitarget  tracking  becomes  clear  with  use. 

The  PMHT  algorithm  assumes  independent  targets.  In  the  special  case  of  linear-Gaussian 
statistics,  each  target  has  a  state  motion  model  of  the  form 


where  w,  is  white  Gaussian  process  noise  with  known  covariance  matrix  Ot-  Measurements  z, 
within  each  scan  Z,  are  assumed  conditionally  independent,  when  conditioned  on  the  collection  of 
target  states.  Measurements  in  different  scans  are  also  assumed  conditionally  independent. 
Different  scans  may  have  different  numbers  of  measurements. 

The  PMHT  observer  is  defined  by  the  state  of  the  postprocessor.  Since  measurements  are 
not  presented  to  the  postprocessor  with  labels  identifying  the  correct  measurement-to-track 
assignments,  the  appropriate  observer  state,  denoted  Ot ,  for  measurement  scan  Z,  comprises  the 
collection  of  target  states  together  with  the  collection  of  measurement-to-track  assignments. 
Therefore,  the  observer  state  has  both  continuous  and  discrete  components  that  must  be 
estimated.  See  the  last  paragraph  in  this  section  for  a  brief  description  of  the  PMHT  estimation 
methodology  (EM). 

Suppose,  momentarily,  that  the  observer  state  is  defined  so  that  it  does  not  comprise 
measurement  assignments.  Then,  even  if  the  target  states  are  known  exactly,  this  modified 
observer  is  unable  to  determine  the  correct  measurement-to-track  assignments.  Hence,  this 
modified  observer  cannot  be  an  observer  in  the  strict  sense  (see  reference  4,  section  9.2). 
Consequently,  the  measurement-to-track  assignments  must  be  part  of  the  observer  state  definition. 
Incorporating  assignments  explicitly  into  the  observer  state  appears  to  be  novel  to  this  study. 

The  probability  density  function  (PDF)  of  the  measurement  scan  Z,  is  conditioned  solely  on 
the  scan  observer  state  Ot.  The  PDF  of  a  measurement  z,  in  Z/,  conditioned  on  Ot,  is  assumed 
known.  In  the  linear-Gaussian  case,  the  measurement  PDF  has  the  form 

where  A/ is  the  discrete  component  of  Ot  corresponding  to  z,;  i.e.,  kt  denotes  the  target  of  origin 
for  z,.  This  is  equivalent  to  writing 


where  v  is  white  Gaussian  measurement  noise  with  known  covariance  matrix  R  .  Because  the 

scan  measurements  are  assumed  to  be  independent,  conditioned  on  the  observer  state,  the  scan 
likelihood  function,  denoted  by  P(Z,|(9,),  is  the  product  over  all  measurements  z,  in  Z,  of  the 
individual  measurement  likelihood  functions  p(z,\Oi).  It  is  important  to  emphasize  that  the  discrete 
components  of  the  observer  O,  are  essential  to  the  definition  of  the  scan-likelihood  function. 

PMHT  is  a  batch  algorithm  that  uses  a  finite  number  of  successive  measurement  scans  to 
estimate  the  batch  observer  state  O.  The  batch  observer  is  the  collection  of  scan  observers  in  the 
batch,  and  the  batch  observer  state  comprises  the  collection  of  the  scan  observer  states,  so  that 
0=  {Ot}.  A  batch  PDF  is  formulated  for  the  batch  observer  using  the  target  motion  models  and 
the  scan  likelihood  functions.  The  batch  PDF  is  a  joint  function  of  the  measurements,  the 
continuous  component  X,  and  the  discrete  component  K  of  the  batch  observer  state  O. 

BINs  are  used  in  this  report  to  display  the  conditional  independence  assumptions  between 
measurements  and  the  batch  observer  states  in  a  graphical  manner.  These  conditional 
independence  assumptions  are  fundamental  to  the  batch  observer  Joint  PDF  and  are  easily 
understood  if  presented  in  a  BIN  graphical  form.  The  BIN  for  PMHT  is  shown  in  figure  1 .  The 
nodes  and  edges  of  the  BIN  graph  are  interpreted  in  such  a  way  that  the  graph  is  equivalent  to  a 
factorization  of  the  batch  observer  PDF.  (See  the  end  of  section  3  for  a  discussion  of  BINs.) 


TIME  t-1 


TIMEt 


Figure  1.  BIN  Fragment  for  PMHT  at  Scan  t 
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The  MAP  estimate  of  O  is  the  maximum  (over  the  observer’s  continuous  and  discrete 
components)  of  the  PDF  of  O  conditioned  on  the  available  batch  measurements.  However,  in  the 
absence  of  computationally  efficient  search  techniques  (e.g.,  dynamic  programming),  the  MAP 
estimate  is  combinatorially  hard  to  compute  because  it  involves  enumeration  of  all  possible  hard 
measurement-to-track  assignments.  Computation  of  such  assignments  is  not  required  by  PMHT. 

The  PMHT  algorithm  differs  from  the  MAP  point  estimate  just  described.  The  PMHT 
estimates  of  the  batch  observer  state  are  computed  in  three  separate,  but  interrelated,  steps.  In 
the  first  step,  the  batch  joint  PDF  is  marginalized  (summed)  over  the  observer  discrete  component 
K,  i.e.,  over  all  measurement-to-track  assignments.  The  marginal  density  is  the  joint  PDF  of  the 
batch  measurement  Z  and  the  continuous  component  X  of  the  batch  observer.  The  second  step 
estimates  the  target  states  x  from  the  marginal  PDF  using  an  algorithm  derived  by  the  EM 
method.  This  computation  yields  the  PMHT  state  estimate  for  each  target  at  each  scan  in  the 
batch.  In  the  third  step,  the  conditional  density  for  the  observer  discrete  component  K  is 
computed  from  Bayes  theorem  by  conditioning  on  the  measurements  Z  and  on  the  PMHT  target 
state  estimates  computed  in  the  second  step.  This  conditional  PDF  is  the  probability  of 
measurement-to-track  assignment  for  every  measurement  and  track  pair. 

The  coupling  between  the  three  PMHT  steps  is  of  fundamental  importance.  The  coupling  is 
caused  by  assignment  interference  because  it  must  estimate  the  parameters  defining  the  discrete 
conditional  distribution  of  the  observer  discrete  component  K.  These  parameters  are  merely  the 
measurement-to-track  assignment  probabilities.  The  PMHT  algorithm  is  essentially  a  recursive 
algorithm  for  estimating  assignment  interference. 

The  second  PMHT  step  yields,  upon  convergence,  a  conditional  error  covariance  matrix  for 
each  target.  It  will  be  shown  that  the  inverse  covariance  matrix  for  target  s,  computed  in  the 
second  step,  is  the  expected  FIM,  where  the  expectation  is,  with  respect  to  the  conditional 
density,  on  the  observer  discrete  component  K. 

PMHT  reliance  on  marginalization  is  unusual  from  the  perspective  of  multitarget  tracking, 
but  it  is  natural  from  a  probabilistic  perspective.  Marginalization  is  also  a  traditional  method  of 
treating  so-called  nuisance  variables  in  statistical  problems.  Marginalization  removes  dependence 
on  knowledge  of  particular  outcomes  of  the  variables  marginalized,  but  the  dependence  on  the 
distributional  parameters  of  these  variables  remains.  In  the  multitarget  tracking  problem,  the 
measurement-to-target  assignments  are  the  nuisance  variables,  and  they  are  marginalized  out. 

The  parameters  of  the  nuisance  variable  distributions  are  estimated  from  the  marginal  distribution 
using  the  EM  method. 
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3.  PMHT  OBSERVER  LIKELIHOOD  STRUCTURE 


Let  A/  >  1  denote  the  assumed  number  of  independent  target  motion  models.  The  integer 
M  should  be  at  least  as  large  as  the  number  of  targets  present  in  the  measurement  sequence. 
Choosing  M  greater  than  the  number  of  targets  may  have  practical  utility  for  background  noise  or 
clutter  normalization  purposes. 

The  PMHT  algorithm  is  a  batch  algorithm.  Let  T>\  denote  the  length  of  the  current 
batch.  The  measurement  scans  are  numbered  so  that  the  batch  comprises  scans  from  time  t  =  \ 
to  time  /  =  T.  The  PMHT  observer  is  a  batch  observer,  so  it  is  convenient  to  define  the  batch 
observer  as  a  collection  of  scan  observers.  For  /  >  1 ,  the  PMHT  scan  observer  state  is  denoted  by 
=  (A'j,  a:,),  and  it  comprises  the  continuous  component  vector  X/  and  the  discrete  component 

vector  Kt-  The  continuous  component  Xt  comprises  the  M  target  states.  Explicitly,  Xt  is  given  by 

where  denotes  the  state  of  target  s  at  time  t.  Component  states  Xo  and  x,m  o^Xt  are 
independent,  for  s^m  because  the  targets  are  assumed  to  be  independent.  The  a  priori  PMHT 
scan  observer  state  is  denoted  by  0^  =  where  the  discrete  component  =0  (the 

empty  set)  because,  by  definition,  measurements  are  unavailable  at  /  =  0.  The  continuous 
component  Afg  =(xg|,...,Xgj^^^  is  discussed  below.  The  target  state  vector  X  comprises  the 
state  of  target  m  at  each  time  in  the  batch,  that  is. 

Let  denote  the  dimension  ofthe  target  state  for  model  5.  The  state  dimension  may  vary 
from  target  model  to  target  model,  but  it  is  assumed  constant  from  scan  to  scan. 

Let  x^  j  denote  the  process  model  for  target  s,  \<s<  M.  Target  model  5  assumes 
that  Xg^  is  a  realization  ofthe  ci priori  distribution  denoted  by  time  dependence  of 

the  functions  j  j  is  indicated  implicitly  by  the  function  argument.  In  the  special  case  of 
linear-Gaussian  target  process  models. 


where  S)  denotes  a  multivariate  Gaussian  PDF  with  mean  vector  p  and  covariance  matrix 
Z.  Equation  (1)  is  equivalent  to  the  conventional  form 


=  F  ,  X  ,  +G  ,  w 

t~\,s  t-\,s  t-],S  t- 


(2) 


where  j  ^  is  a  white  Gaussian  process  with  (positive  definite)  covariance  matrix  0^  ^  ^ .  Noise 
processes  and  n',.,  ^  are  assumed  to  be  independent  because  5  w.  The  matrices  F  ^  ^ , 

G  ,  ,  and  0  ,  .are  assumed  known  for  all  /  >  0  and  s. 

/-1, 5’  -"/-1, 5’ 


The  target  process  models  hold  for  /  =  0,  so  state  is  a  realization  of  the  distribution 


^X^uK)  ^s(^0s)^0s 


(3) 


The  initial  target  states  can  therefore  be  eliminated  without  altering  the  batch  likelihood 

structure.  In  this  report,  however,  target  states  are  initialized  at  time  /  =  0  for  notational 
convenience.  For  the  linear-Gaussian  case,  x^^^  is  a  realization  of  the  a  priori  Gaussian  PDF 


(4) 


where  the  mean  x^^  and  covariance  matrix  are  given,  and  the  distribution  (3)  is  given  by 
(p  fx  )=  A^(x,  F  x„  ,  F  F'  +G„  0,  G'  ), 

'  \s J  \  U  Oj  0^’  Oj  ^0^  Oi-  os/^ 

where  the  mean  and  covariance  of  J  are  the  predictive  statistics. 


Let  denote  the  PDF  of  the  scan  observer  continuous  component  X,.  The  time 
dependence  of  will  be  indicated  implicitly  by  its  argument.  For  /  >  1,  the  PDF  of  A)  is 
conditioned  on  X  .  From  the  assumption  of  independent  targets,  it  follows  that 


M 


(5) 


The  a  priori  PDF  of  the  continuous  component  Xo  defining  the  initial  observer  is  given  by 

'i'.w=n 

(6) 


where  independent  target  initialization  has  been  assumed. 
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For  ?  >  1 ,  the  discrete  component  Kt  of  the  scan  observer  state  O,  comprises  a  unique 
measurement-to-track  assignment  for  every  measurement  in  scan  Z,.  The  discrete  component  Kt 
and  the  continuous  component  Xt  are  assumed  to  be  statistically  independent.  Let  a?,  denote  the 
number  of  measurements  in  scan  Z,.  For  notational  simplicity  it  is  assumed  that  a?,  >  1 ;  however, 
no  theoretical  difficulty  arises  if  aa,  =  0.  Recall  that  Ko  =  0.  Let  Ztr  denote  measurement  r  in  scan 

Z^  =  (z^, ,  , . . . , Zj  ^  j .  The  discrete  component  Kt  is  defined  by 


where  1  <  ^  A/  for  all  a*  =  1 ,2, . . .  ,a/,.  Thus,  measurement  Ztr  is  assigned  by  Kt  to  the  track  motion 

model  with  index  eK^.  It  is  assumed  that  the  discrete  components  |  are  independent,  that 

is,  A, rand  are  independent  if  r  r'.  Finally,  the  discrete  component  vectors  J  are  assumed 

statistically  independent  from  scan  to  scan,  that  is,  K  is  independent  of  ATr, because  t^t'. 

The  PMHT  measurement  PDF  is  defined  for  measurement  z,r  by 

All  measurements  in  all  scans  are  assumed  (for  simplicity)  to  have  the  same  dimension,  denoted  by 
Nz-  Under  linear-Gaussian  assumptions, 

C  (z  lx  )  =  ‘^(z  \H  X  ,R  ).  (7) 

^  m\  tr\  tm  /  \  tr\  tm  tm  /  ^ 

The  linear-Gaussian  case  can  also  be  written  in  the  equivalent  form 

^tr=».m^tn.+^t.^ 

where  v^,,  is  an  additive  white  Gaussian  noise  process  with  covariance  Rtm-  The  PMHT  scan 
likelihood  function  is  a  PDF  defined  over  the  measurements  in  scan  Z,,  and  it  is  conditioned  on  the 
observer  state  Ot.  Explicitly,  because  the  measurements  are  conditionally  independent, 

p{z\p).p{z\X„K).Y{ 

-1  '  '  (8) 

The  discrete  component  k^  of  the  vector  AT,  serves  as  a  pointer  to  the  appropriate  target  model 
and  is  thus  an  integral  component  of  measurement  conditioning.  The  parameters  defining 
probability  mass  function  of  the  discrete  components  of  the  measurements  in  scan  Kt  are  not  part 
of  the  observer  state.  (Reference  to  the  BIN  of  figure  1  will  clarify  this  important  distinction.) 

Let  represent  the  probability  that  a  measurement  in  scan  Z,  is  assigned  to  target  motion 
model  AAA,  where  /  >  1 .  The  probability  reflects  the  fraction  of  scan  measurements  assigned  to 
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target  m  at  time  t.  These  target  measurement  probabilities  are  needed  because  some  targets  may 
produce  more  measurements  per  scan  than  others  because  of  individual  target  characteristics 
(e.g.,  signal-to-noise  ratio),  environmental  effects,  sensor  properties,  and  other  application 
considerations.  Denote  the  "within-scan"  measurement  probability  vector  by 

(9) 

The  vector  ;r,  parameterizes  the  distribution  of  the  discrete  component  k,r  for  the  measurement  z,r 

in  scan  Z,.  Because  measurements  within  a  scan  are  assumed  conditionally  independent  and 
identically  distributed,  the  distribution  ;r,  is  assumed  to  be  the  same  for  all  measurements  made  at 

time  t.  Thus,  =  Probjit,^  =  wj  for  all  measurement  assignments  k,r.  The  batch  target 

measurement  probabilities  IT  =  11^.  =  are  estimated  by  the  PMHT  algorithm  from  the 

batch  measurement  data..  Alternative  assumptions  concerning  IT  are  possible.  For  instance,  a 
Bayesian  a  priori  density  can  be  assumed  for  if  desired.  A  Bayesian  a  priori  density  for  Kt  is 

analogous  to  the  initialization  (3)  for  A,.  This  and  other  alternatives  are  described  in  section  7. 1 . 

Let  denote  the  discrete  PDF  (or,  probability  mass  function)  of  the  discrete  component 
Kt  of  the  scan  observer  Ot.  The  time  dependence  of  'Fjj.  is  indicated  implicitly  by  its  argument. 
The  PDF  of  K,  takes  the  form 


'r,(^,)=Prob 


(10) 


The  product  in  (10)  follows  from  the  assumption  that  different  measurements  within  a  scan  have 
statistically  independent  and  identically  distributed  discrete  components  ktr.  It  turns  out  (see 
equation  (22)  and  the  ensuing  discussion)  that  the  measurement  independence  assumptions  imply 
that  Ktm  are  mixing  proportions  of  a  mixture  density  modeling  the  scan  measurements.  For  the 
linear-Gaussian  case,  the  means  of  the  Gaussian  components  of  the  mixture  are  the  target  state 
estimates  for  the  scan. 


The  batch  observer  state,  denoted  0  =  0^=  comprises  the  continuous  and 

discrete  components  of  the  T-scan  observers  in  the  current  batch.  The  dependence  of  the 
observer  and  its  components  on  batch  length  Tis  suppressed  throughout  the  remainder  of  this 
report.  Thus,  the  continuous  and  discrete  components  of  the  batch  observer  O  are 


and 


(11) 


(12) 
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respectively.  is  a  place  holder  in  (12)  because  =  0-  The  batch  measurement  is  denoted 


It  is  assumed  that  measurement  scans  and  measurements  within  scans  are  statistically  independent, 
conditioned  on  the  batch  observer  state  O. 

The  PMHT  batch  observer  PDF  is  a  joint  function  of  the  measurements  Z  and  the  batch 
observer  state  O.  Because  no  measurements  are  made  at  time  /  =  0,  the  a  priori  PDF  of  the  scan 
observer  at  time  /  =  0  is  defined  to  be 

A/ 

'*'(o.)=n?’.k.)  <1“) 


The  scan  observer  at  time  /  >  1  is  conditionally  dependent  on  the  scan  observer  at  time  /  -  1 .  The 
scan  observer  continuous  and  discrete  components  are  independent,  so  the  observer  conditional 
PDF  is  given  by 


Substituting  the  expressions  (5)  and  (10)  gives  for  t  >  1 


The  batch  observer  PDF  is,  from  the  conditional  independence  assumptions. 


(15) 


Piz,o)  s  Piz.K  K) = >F(o.)n  'p(o,|o,.,)f(z,|o,), 

f=l 


(16) 


Substituting  (8)  and  (15)  into  this  expression  gives  the  batch  observer  joint  PDF  in  the  form 
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When  the  parametric  dependence  on  the  probabilities  n  is  made  explicit,  the  joint  PDF  is  written 
P{Z,X,K :  n).  The  batch  PDF  (17)  is  the  fundamental  probabilistic  structure  of  the  PMHT 
algorithm.  The  derivation  of  the  PMHT  algorithm  is  given  in  section  4. 

Insight  and  understanding  of  the  structure  of  P{Z,X,K)  is  enhanced  by  displaying  its 
conditional  independence  assumptions  in  a  mathematically  equivalent  graphical  format.  A  BIN  is 
a  very  effective  representational  technique  designed  for  just  such  a  purpose.  A  fragment  of  the 
fundamental  graphical  structure  of  the  PMHT  batch  PDF  for  measurements  at  scan  t  is  portrayed 
as  a  BIN  in  figure  1 .  The  BIN  of  figure  1  is  a  directed  graph,  and  its  nodes  represent  random 
variables  in  the  batch  observer  PDF  P{Z,X,K).  The  directed  edges  depict  the  conditioning  of 
these  random  variables  in  the  following  way:  Each  node  is  conditioned  jointly  on  all  of  its 
"parents. "  The  overall  likelihood  structure  of  the  BIN  is  defined  to  be  the  product  of  all  the 
conditional  PDFs.  There  are  precisely  as  many  conditional  PDF  factors  as  there  are  nodes  in  the 
graph.  Following  this  procedure  for  the  BIN  of  figure  1  shows  that  its  joint  likelihood  function  is 
identical,  factor  by  factor,  to  the  factorization  of  P{Z,X,K)  given  by  equation  (17).  The  precise 
mathematical  forms  of  the  individual  conditional  PDF  factors  are  not  defined  by  the  BIN,  but  they 
are  stipulated  by  the  requirements  of  the  specific  application.  The  Gauss-Markov  equation  (2), 
the  Gaussian  initialization  (3),  and  the  measurement  conditioning  equation  (7)  are  specific  to 
PMHT.  The  directed  graph  provides  global  likelihood  structure,  while  the  specific  conditional 
PDFs  provide  local  structure.  The  graph  may  not  have  "directed  cycles"  because  they  result  in 
improper  factorizations  of  the  global  joint  likelihood  function. 

The  special  caseM=  1  of  the  PMHT  batch  PDF  (17)  is  identical  to  the  PDF  of  a  fixed- 
interval  Kalman  filter  formulated  for  multiple  measurements.  The  reason  for  these  identical  PDFs 
is  that,  in  this  case,  the  batch  discrete  component  K  and  the  within-scan  target  measurement 
probabilities  H  are  trivial  because  there  is  only  one  target  and  all  batch  measurements  are 
necessarily  assigned  to  it.  In  addition,  if  each  scan  has  precisely  one  measurement,  the  PMHT 
batch  likelihood  function  is  identical  to  the  likelihood  function  of  the  usual  textbook  Kalman  filter 
(see  reference  5). 
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4.  DERIVATION  OF  THE  PMHT  ALGORITHM 


In  this  discussion,  the  authors  assume  the  reader  is  familiar  with  the  EM  method.  For 
background  material  on  this  very  general  technique,  see  reference  6  and  the  references  mentioned 
therein.  For  PMHT,  the  "missing  data"  in  the  sense  of  EM  is  the  discrete  component  K  of  the 
batch  observer,  and  the  “complete  data”  comprises  Z,  X,  and  K.  This  confusing  nomenclature  is 
not  used  in  this  report.  In  the  special  case  of  linear-Gaussian  statistics  with  known  covariances, 
the  target  measurement  probabilities  IT  and  the  states  X  of  the  continuous  component  are  batch 
parameters  to  be  estimated.  This  case  is  especially  simple  because  the  parameters  defining  the 
target  process  and  measurement  PDFs  are  linear  fonctions  of  the  state  means,  and  the  estimated 
means  are  the  target  state  estimates.  In  general,  however,  it  is  necessary  to  distinguish  between 
state  parameter  estimates  and  state  estimates.  For  nonlinear,  non-Gaussian  processes,  therefore, 
it  is  assumed  that  X  is  the  parameter  set  defining  the  target  process  and  measurement  PDF.  This 
slight  abuse  of  notation  should  not  cause  confusion. 

The  discussion  in  this  section  derives  the  expectation  step  (E-step)  of  the  EM  method  for  the 
general  case  because  the  E-step  is  done  as  easily  in  general  as  it  is  for  linear-Gaussian 
assumptions.  As  is  seen  below,  the  measurement  probabilities  11  and  target  parameters  X  are 
treated  separately  in  the  maximization  step  (M-step)  even  in  the  general  case.  The  M-step  treats 
the  probabilities  n  for  the  general  case.  The  M-step  for  the  parameter  estimates  X  is  significantly 
easier  to  understand  for  linear-Gaussian  statistics,  so  this  special  case  is  treated  separately  first. 
Subsequently,  the  general  case  for  X  is  treated. 

The  E-step  begins  by  defining  a  PDF  on  the  discrete  component  K.  Let  denote  this  PDF. 
On  convergence  of  the  PMHT  algorithm,  ^  gives  the  measurement-to-track  assignment 
probabilities  for  all  possible  measurement  and  target  pairs.  From  Bayes  theorem,  the  conditional 
PDF  on  K  is  defined  by 


^(k\z,x) 


p{z,x,  K  n) 
P(Z,  A :  n)  ’ 


(18) 


where  the  denominator  of  (18)  is  the  marginal  distribution  of  the  PMHT  likelihood  function  over 
the  discrete  component  K.  When  the  parametric  dependence  on  the  probabilities  IT  is  made 
explicit,  the  conditional  PDF  (18)  is  written  x{K\Z,Xlt^  The  marginal  distribution  of  P(Z,X,K) 
over  K  is  defined  by 


P(Z,A)  =  P(Z,A:n)s2;  P^2,X,K), 

K 


(19) 


where  the  summation  over  the  batch  component  K  is  defined  by 

H-tY.  -Z-S. 

K  K,  K,  Kj. 


(20) 
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and  where  the  sum  over  the  scan  component  Kt  is  defined  by 


"/  M  M  M 

S  I  -z  -z . 


(21) 


Using  the  most  expanded  of  the  summation  forms  (20)  and  (21),  it  is  straightforward  to  verify  the 
important  algebraic  identity 


P{Z,X)  = 


1 

1 

n 

r=l 

V  ^  u  X  ) 

Z— im  ^  m\  tr\  tm/ 
_m-\ 

(22) 


PMHT  computes  an  ML  estimate  of  the  continuous  component  X  and  the  target  measurement 
probabilities  IT  from  the  marginal  PDF  (22).  The  distribution  P{Z,X)  is,  thus,  the  appropriate 
PDF  for  a  batch  observer  whose  state  definition  does  not  explicitly  include  discrete  components. 
The  marginal  PDF  is  fundamental  to  multitarget  observability  questions,  and  observability  is 
related  to  PMHT  algorithm  convergence.  (See  section  6  for  additional  information.) 

The  marginal  PDF  (22)  exhibits  clearly  the  underlying  assumption  that  the  measurements 
within  scan  are  conditionally  independent.  It  also  provides  the  useful  interpretation  that  the 

components  of  the  batch  target  measurement  probability  vector  n,  are  the  mixing  proportions  of 
the  measurement  PDF  mixture 


from  which  the  measurements  in  scan  Z,  are  drawn.  Substituting  (17)  and  (22)  into  definition  (18) 
gives 


X(,K\Z,X)  =  f[  fi 


w 


(23) 


where  the  weight  >  0  is  a  function  of  Z,  X,  and  fl  and  is  given  by 


vv  = 

str 


TT  C  \x  ] 

Is^  ^  \  t/'  1  ts  / 

iz,  X  ) 

Lma  tm~  m  \  /r  1  tm  / 

m=] 


(24) 


14 


The  expression  (23)  does  not  include  a  product  over  the  state  models  because  terms  involving 
cancel  out  of  the  ratio  (18).  The  weight  is  interpreted,  using  Bayes  theorem,  as  the 

conditional  probability  that  measurement  Ztr  is  assigned  to  the  target  model  s,  conditioned  on  the 
continuous  component  X  and  the  measurements  Z;  that  is, 

w,^  =  Prob[*,^=5|Z,x]. 


From  the  definition  (24),  it  is  straight  forward  to  verify  the  algebraic  identities 
2;  :^(K\Z,X)=l, 

K 

and 


(25) 


X  X{K\Z,X)^w^^^ 


(26) 


where,  in  (26),  the  sum  over  K\k,r  is  the  sum  over  all  indices  in  K  except  K. 

The  E-step  is  concluded  by  evaluating  the  expectation  of  the  logarithm  of  the  PMHT  batch 
PDF  (17),  where  the  expectation  is  with  respect  to  the  conditional  PDF  (23).  Let  X'  denote  a 
given  value  for  the  continuous  component  of  the  batch  observer,  and  let  11'  denote  a  given  value 
for  the  target  measurement  probabilities.  The  primes  on  these  variables  do  not  denote  the 
matrix/vector  transpose  operator.  The  expectation  required  by  the  E-step  is  written  explicitly  as 


Q^Q(U,X\U’,X')  =  Y.  {\ogP{Z,X,K.Yl)]  7:{K\Z,X':Yl'). 

K 


(27) 


The  function  Q  is  called  the  auxiliary  function  in  reference  6,  and  it  is  closely  related  to  cross¬ 
entropy  and  the  Kullback-Leibler  distance.  In  this  report  function  Q  is  called  the  cross-entropy 
function.  Taking  the  logarithm  of  the  PMHT  batch  PDF  (17)  gives 

log  P(Z,X,K:n)  =  I;  log  +  IZ  los  9,KKu) 

v=i  r=i  j=l 

+  +  log 

^=1  r~] 

Substituting  this  expression  into  the  definition  (27),  interchanging  the  summation  order  so  that, 
for  example, 

z  z  z  =z  z  z  z . 

Kir  >  ' 
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and  then  using  the  identities  (25)  and  (26)  gives  the  cross-entropy  function  in  the  form 

T  M 

e  =  I  e,.n  + 1 

r=l  m=l 

where 

a,n  •  =  i  t 

r=\  m=\ 

Q..X  • 

=  log  (i’.KJ+zjiog  •p.(’‘,.K,J+'Z  log  f.(z„k„)|, 

and  where  the  weights  =  w'^^^[Z,X',n')  are  defined  as  in  (24).  In  equations  (29)  and  (30), 
the  notations 

and 

have  been  used.  Again,  the  primes  on  these  variables  do  not  denote  transpose. 

The  M-step  is  a  maximization  problem.  Explicitly,  given  values  O'  and  the  M-step 
requires  computing  values  n  and  X  for  which 

e(n,i'|n',j^')  =  max  Q(u,x\n’,x'). 

(31) 

The  given  variables  O'  and  X'  comprise  the  current  values  of  the  PMHT  algorithm,  and  the 
variables  n  and  X  comprise  updated  values  of  the  algorithm.  The  PMHT  recursion  is  stated 
explicitly  in  section  5.  On  convergence  of  the  PMHT  algorithm.  Hand  X  are  ML-parameter 
estimates.  Equations  (28),  (29),  and  (30)  demonstrate  that  the  maximization  problem  (31) 
decouples  into  a  maximization  problem  for  each  of  the  /’-probability  vectors  tt,  and  a 
maximization  problem  for  each  oftheM-target  state  sequences  X" .  The  weights  change 
from  M-step  to  M-step,  but  within  each  M-step,  the  weights  are  fixed. 


(28) 


(29) 


(30) 
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The  maximization  problem  for  itt  is  constrained  by  the  requirement 


M 


Z 


n. 


=  1, 


(32) 


so  the  appropriate  Lagrangian  for  this  problem  is 


M 


m-l 


where  yt  is  the  Lagrange  multiplier.  Differentiating  the  Lagrangian  with  respect  to  rttm  and  setting 
the  result  to  zero  gives  the  necessary  condition 


tm  LmU  mtr 

r,  r=l 


Summing  these  equations  from m  =  I  toM and  using  the  constraint  (32)  gives 


M 


r.=I.  Z 

r=]  m=l 


W 


From  the  expression  (24)  for  the  weights,  it  follows  easily  that  y  =  n,,  so  the  unique  stationary 
point  of  the  Lagrangian  is 


w'  . 

mtr 


(33) 


By  lemma  2  of  reference  7,  it  follows  that  ^  =  is  the  unique  global  maximum  of 

the  cross-entropy  function  It  follows  from  (33)  and  (24)  that  =  0 ,  if  and  only  if 
k,  =  0 .  Because  Q  n  =  and,  hence,  O  =  -oo  if  ^  =  0 ,  it  is  assumed  without  loss  of 
generality  that  the  initial  probabilities  H'  are  chosen  to  be  strictly  positive. 

The  solution  for  the  state  sequence  X”  for  target  m  is  derived  first  for  the  special  case  of 
linear-Gaussian  statistics.  The  insight  provided  by  this  special  case  is  helpful  in  understanding  the 
general  case,  which  is  discussed  beginning  with  equation  (39).  Recall  the  general  gradient  identity 


(Fx-^yi-^Fx-M)  =  2F’l-\Fx-n). 


(34) 
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Substituting  the  linear-Gaussian  models  (1),  (4),  and  (7)  into  the  definition  of  the  cross-entropy 
fianction  and  taking  the  gradient  of  Q„,x  with  respect  to  the  state  vector  x,™  for  /  = 
gives  a  symmetric  block  tridiagonal  system  of  linear  equations  for  the  state  sequence  for  target 
model  m.  This  system  is  written 


(35) 


where  the  block  matrices  Ap„,  Dtm,  and  Btm  are  given  by 
A  y~' 

Om  ^Om 

B  =F'(G  Q  G'  V,  0<t<T-\ 

tm  tm\  tm^tm  tm)  ^  ^  a, 

and  the  synthetic  measurement  is  defined  by 


nn, 

/  tm 


z  , 

mtr  tr  ’ 


r=\ 


\<t<T. 


(36) 


(37) 


The  solution  of  the  system  (35)  is  the  updated  M-step  state  sequence  X” .  This  completes  the 
algorithm  derivation  for  the  linear-Gaussian  case. 

The  number  in  (36)  and  (37)  represents  the  expected  number  of  measurements  in  scan 
Z,  that  are  assigned  to  target  m.  From  (33),  it  follows  that  the  synthetic  measurement  is  the 
probabilistic  centroid  for  target  m  of  the  measurements  in  scan  Z,;  that  is,  zj^  is  the  expected 
measurement  for  target  m  at  time  t.  The  probabilistic  centroid  z^^  always  lies  in  the  convex  hull 

of  the  scan  measurements.  (The  convex  hull  of  a  given  set  is,  intuitively,  the  smallest  convex  set 
that  contains  the  given  set  as  a  subset.  By  definition,  the  convex  hull  is  the  intersection  of  all 
convex  sets  containing  the  given  set.) 

The  solution  of  the  block  tridiagonal  system  (35)  can  proceed  along  strictly  algebraic  lines, 
following  the  methods  suggested  in  section  5.5  of  reference  8.  Unfortunately,  this  procedure 
provides  little  insight  into  the  resulting  algorithm,  and  it  does  not  yield  state  error-covariance 
matrices.  A  direct  connection  with  available  Kalman  filtering  techniques  overcomes  both  these 
deficiencies. 
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The  connection  is  established  by  noting  that 
Kalman  filter.  From  (30), 


is  the  PDF  of  a  linear-Gaussian 


(38) 


where  the  target  measurement  centroid  is  defined  by  (37).  These  proportional  expressions  are 

derived  by  algebraically  manipulating  the  state-dependent  terms  of  the  Gaussian  exponents.  The 
right-hand  side  of  (38)  is  the  PDF  of  a  Kalman  filter  whose  plant  is  identical  to  the  target  m 
process  model  and  whose  measurement  is  the  target  measurement  centroid  .  The  conditional 

independence  relationships  implicit  in  likelihood  function  (38)  are  displayed  graphically  in  the  BIN 
of  figure  2.  The  directed  graph  of  figure  2  is  a  subgraph  of  the  full  PMHT  graph,  where  the  full 
graph  comprises  all  the  BIN  fragments  depicted  in  figure  1 .  Well-known,  fixed-interval  Kalman 
filter-smoothing  recursions  (reference  5,  section  7.4)  can  therefore  be  used  to  compute  the  ML- 
state  sequence.  These  recursions  are  used  in  the  PMHT  algorithm  summary  described  below  in 
equations  (49)  through  (56).  This  approach  shows  that  the  solution  of  the  block  tridiagonal 
system  (35)  is  identical  to  the  state  estimates  of  a  fixed-interval  Kalman-smoothing  filter. 

If  the  solution  of  (35)  or  (37)  exists  and  is  unique,  then  it  maximizes  the  cross-entropy 
function  as  required  by  the  M-step.  The  fact  that  such  a  solution  must  be  a  maximum  for 

is  not  evident  from  the  strictly  algebraic  development,  but  it  follows  immediately  from 

equation  PDF  (38).  In  general,  however,  the  solution  may  not  exist  and,  if  it  exists,  it  may  not  be 
unique.  Significant  insight  into  these  possibilities  is  provided  by  the  PDF  (38)  because  they  are 
related  to  the  observability  problems  associated  with  the  Kalman  filter.  Further  discussion  of 
convergence  issues  are  given  in  section  6. 


Figure  2.  BIN  Structure  of  an  Interval  Kalman  Filter 
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Error  covariances  corresponding  to  the  state  sequence  X""  are  readily  computed  using 
established  fixed-interval  recursions  for  the  Kalman-smoothing  filter  (38).  These  recursions  are 
given  by  (57)  and  (58).  For  further  discussion  and  derivation  of  these  recursions,  see  page  189  of 
reference  5  and  the  references  therein.  It  is  stressed,  however,  that  the  M-step  of  the  PMHT 
algorithm  does  not  explicitly  use  error  covariances. 


It  is  intuitively  clear  that  the  PMHT  error  covariances  computed  from  (38)  for  X”  are 
important  in  the  multitarget  tracking  application.  However,  the  cross-entropy  function  does  not 
directly  provide  a  statistical  interpretation  for  error  covariances.  Section  7  gives  such  an 
interpretation  in  terms  of  a  randomized  decision  rule  applied  to  the  measurement-to-track 
assignments.  As  will  be  seen,  this  interpretation  also  holds  in  the  nonlinear  non-Gaussian  case. 


The  M-step  for  the  PMHT  algorithm  is  now  easily  treated  for  the  general  case.  As  discussed 
in  the  beginning  of  this  section,  for  nonlinear  non-Gaussian  processes,  X  denotes  a  parameter  set 
defining  target  and  measurement  distributions,  and  not  state  estimates  as  in  the  linear-Gaussian 
case.  Instead  of  (38),  one  obtains  from  (30)  the  expression 


X 


(39) 


(40) 


where,  for  an  appropriate  normalization  constant  c. 


c  fl  icjzjxjr" 


(41) 


defines  a  conditional  PDF  on  the  full-scan  measurement.  (The  conditional  PDF  (41)  is  different 
from  the  observer  scan  PDF  (8).)  The  parameter  sequence  X"  that  maximizes  (40)  is  the 

updated  parameter  set  required  by  the  M-step.  The  computation  of  A”  requires,  in  general, 
using  an  iterative  numerical  algorithm.  This  numerical  procedure  is  conceptually  equivalent  to  a 
single-target  MAP  tracker,  and  its  availability  is  assumed.  Whether  the  MAP  tracker  is 
computationally  efficient  is  irrelevant  to  the  M-step:  All  that  is  necessary  is  that  it  compute  the 

MAP  parameter  estimated” .  The  computation  of  A"  using  the  MAP  tracker  completes  the  M- 
step  in  the  general  case. 


The  conditional  PDF  (41)  is  assumed  known  because  it  is  either  known  explicitly  or  because 
it  can  be  derived  from  quantities  that  are  given.  In  the  most  general  case,  it  is  possible  to  write 
only 


). 


(42) 
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where  the  measurement  function  htm  is  given,  and  the  PDF  of  the  measurement  noise  v,„  is 
assumed  known.  Theoretically,  the  PDF  can  be  derived  from  (42)  and  the  noise  PDF. 

The  conditional  PDF  j  can  then  be  derived  directly  from  definition  (41). 


Expressions  simpler  than  (41)  are  not  available  in  general.  If  sufficient  statistics  for  the 
samples  are  known  for  the  family  of  distributions  | ,  then  it  is  possible  to  write  the 

product  (41)  succinctly.  A  special  case  is  that  of  additive  Gaussian  measurement  noise,  that  is. 


(43) 


where  denotes  white  Gaussian  noise  with  covariance  Rtm.  Because  (43)  is  equivalent  to 
C  (z  )=  7l\z  h  (x  \R 

~  m  \  tr  /  \  tr  tm\  tm  P  tr 


(44) 


the  PDF  (41)  can  be  written 

E  (Zjx  )  =  7t\l  (Z\\h  (x  ),(nk,  V R^  , 

m\  t\  tm /  \  tm\  t/\  tm\  tm/^\  t  tm)  tm  ^ 

where  the  centroid  is  given  by  equation  (37).  The  sufficient  statistics  are  the 

measurement  centroids  and  the  weighted  covariance  matrices 


The  posterior  state  PDF  corresponding  to  the  parameter  sequence  X"  can  be  derived  from 
the  right-hand  side  of  (40)  and,  subsequently,  state  estimates  derived  from  the  posterior  PDF. 
Error  covariances  may  also  be  derived  from  the  posterior  PDF.  However,  useful  expressions  for 
the  state  estimates  in  the  general  case  are  unavailable. 
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5.  EXPLICIT  STATEMENT  OF  THE  PMHT  ALGORITHM 
IN  RECURSIVE  FORM 


5.1  LINEAR-GAUSSIAN  CASE 


The  PMHT  algorithm  is  summarized  in  this  section  for  the  special  case  of  linear-Gaussian 
statistics.  The  batch  measurement  Z  is  assumed  given.  Initialize  the  target  measurement 

probabilities  so  that  >0.  Initialize  a  target  state  sequence  . 

for  each  of  the  M  target  models.  Let  /  >  0  denote  the  PMHT  iteration  index. 


For  ;■  >  0,  compute  the  assignment  weights 


;r«  \H,  xf\R  ) 

tm  ^  tr '  tm  tm*  Im* 


M 


\H  X^'\R  ) 

ty  '  ts  ts  '  ts> 


s=\ 


(46) 


for  m  =  and  r  =  Thus,  an  assignment  weight  is  computed  for  every 

target  and  measurement  combination  at  each  scan  in  the  batch.  Update  the  target  measurement 
probabilities 


Tt 


0-1) 

tm 


mtr 


(47) 


A  target  measurement  probability  is  computed  for  every  target  and  scan  in  the  batch.  Update  the 
target  measurement  centroids 


5-0*1)  ^ 
tm 


t  tm 


z  , 

mtr  tr  ’ 


t-  \ ...  T 


r=l 


(48) 


A  centroid  is  computed  for  each  target  at  each  scan  in  the  batch.  The  effective  covariance  matrix 
for  the  centroid  (48)  is  proportional  to  the  covariance  matrix  explicitly. 


tm  ^  t  tm  ^  tr, 


t  =  l,...,T. 


(49) 


Target  state  sequences  for  each  of  the  Af  targets  are  updated  using  fixed-interval  Kalman 
smoothers  whose  inputs  are  the  target  measurement  centroids.  Let  the  output  of  the  Kalman 

smoother  for  target  model  m  at  iteration  /  +  1  be  denoted  by  . •  These 

estimates  are  computed  via  a  forward  and  a  backward  recursion.  Initialize  intermediate  variables 
of  the  forward  recursion  by 
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(50) 


(51) 


Intermediate  variables  are  dummy  variables,  so  their  dependence  on  model  m  and  iteration  number 
/■  is  suppressed  for  notational  simplicity.  The  forward  recursion  is  defined  for  t  =  0,1,. ..,r -  1  by 


P  =F  P^F'  +G  Q  G'  , 

t+\\t  tm  //  tm  tm^ 


(52) 


P  zzP  -w  H  P 

/+1'/+1  ;+I|/  l+l  /+l,m  J+l|(> 


y  =F  y  +fV  F  y 

+  l  tm-^  t\t  /+!  /  +  ],m  t+\,m  tnf^  l\l 


(53) 

(54 

(55 


The  updated  PMHT  state  estimate  for  model  m  at  time  T  is  given  by 


f-’)  =  y 

Tm  ^T\T 


(56) 


and  the  updated  PMHT  state  estimates  for  /  =  T-  1,...,1,0  are  given  by  the  backward  recursion 


-  r, 


y  +P  F'  P'^  -F  V  1 

/|(  (m  /  +  lil[  (  +  l,m 


(57) 


Equations  (50)  through  (57)  comprise  a  bank  of  M  Kalman  smoothers  that  are  run  in  parallel  for 
each  PMHT  iteration;  however,  these  smoothing  filters  are  not  independent  because  they  are 

linked  via  the  assignment  weights  defined  by  (46). 


Equations  (46)  through  (57)  comprise  one  step  of  the  PMHT  algorithm  for  the  special  case 
of  linear-Gaussian  statistics.  The  associated  error  covariance  estimates  are  readily  computed. 
Denote  the  PMHT  error  covariance  matrix  of  target  m  at  time  t  by  At  time  T,  =  P^^^. 

The  covariances  at  other  times  in  the  batch  are  computed  by  a  backward  recursion.  Explicitly, 
because  t  =  T-  1,...,0, 


Z  =  P  +  P  F'  P~'  \t  -P 

lit  tit  /m  +  (+11/ 


p-[  F  P, . 

/  +  1|/  im  t\t 


(58) 


The  covariances  P^^  and  used  in  (58)  are  the  intermediate  variables  computed  from  (52)  and 
(54)  during  the  final  PMHT  iteration;  hence,  they  depend  on  the  target  model  and  are  different 
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for  different  targets.  Because  of  the  computational  complexity  of  the  error  covariances  (58),  it  is 
fortunate  that  they  must  be  computed  only  once,  after  the  PMHT  algorithm  has  converged. 

Unlike  the  single-target  Kalman  filter,  the  error  covariances  (58)  cannot  be  computed 
offline.  This  result  is  attributed  to  the  fact  that  Kalman  gains  (53)  are  coupled  by  assignment 
interference;  that  is,  the  number  of  measurements  each  target  receives  is  unknown  and  must  be 
estimated.  The  ML  estimate  of  this  number  is  the  reciprocal  of  the  coefficient  of  Rtm  in  equation 
(49).  It  is  possible  to  compute  the  covariances  offline  only  if  this  coefficient  is  known.  The 
single-target  Kalman  filter  is  thus  a  degenerate  case;  that  is,  the  assumption  of  precisely  one  target 
implies  perfect  knowledge  of  all  measurement  assignments. 

Numerical  problems  can  arise  if  one  or  more  of  the  target  measurement  probabilities 
become  small  during  PMHT  iteration.  This  problem  is  easily  eliminated  by  revwiting  the 

equations  to  cancel  the  common  factors  of  implicit  in  (48).  Explicitly,  the  modified  weights 
are  defined  as 


CO 


(>+i)  ^ , 

mtr  M 


^(2  \H  x^\R  ) 

^  tr'  tm  tm ' 


S=\ 


(59) 


The  updated  target  measurement  probability  ;t  takes  the  form 


(f+l)  — (/-i-l)  (i) 

k\  ^  -  or  ,  V  \ 

tm  mt  tm  ’ 


where  the  mean  measurement  weight  for  target  m  at  time  t  is  defined  by 


mt  mtr 

",  r=l 


The  target  measurement  centroid  defined  by  (48)  is  then  rewritten  in  the  form 


(60) 


(61) 


jo*i)  ^  i —  y 


Similar  modification  of  equation  (52)  gives  the  algebraically  equivalent  form 
fV  =  n  P  (n  ,  P  +P  ,  ) 

f+1  /  +  !  /+l,m  f  +  l|f  f  +  l|m  I  t+\  /  +  r  +  t  +  \,m  /  + 


(62) 


(63) 


The  use  of  equations  (59)  through  (63)  eliminates  numerical  difficulties  associated  with  small 
target  measurement  probabilities. 
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Numerical  difficulties  may  be  encountered  with  generic  Kalman  filtering  if  covariance 
matrices  are  computed  explicitly.  These  kinds  of  difficulties  are  very  closely  related  to  the  well- 
known  numerical  problems  that  can  arise  with  linear  least-squares  computations  if  the  normal 
equations  are  solved  directly.  Strictly  numerical  techniques  such  as  QR  factorization  or  singular 
value  decomposition  (SVD)  (see  reference  8)  are  very  useful  in  practical  computations;  however, 
good  numerical  methods  cannot  prevent  problems  that  arise  because  of  observability  limitations. 
Numerically  stable  forms  of  the  Kalman  filter  have  been  studied.  For  example,  see  reference  9. 


5.2  GENERAL  CASE 


Initialize  the  PMHT  target  measurement  probabilities  11^°^  >0  and  a  target  parameter 
sequence  for  each  of  the  M-targets.  For  /  ^O,  compute  the  assignment 

weights 


(i  +  !) 


4"  (z 

_  ~  m  ^  /r  *  fm  ^ 


(0> 


M 


ts  tr'  ts  ^ 


5=1 


(64) 


for  /M  -  1,...,M,  t  -  l,...,r,  and  r  -  1,...,/?/.  The  target  measurement  probabilities  are  updated 
exactly  as  before  using  equation  (47).  Define  the  conditional  PDF 


r=l 


(65) 


where  c  is  an  irrelevant  normalization  constant.  The  parameters  are  unknown  variables 
in  (65).  Apply  the  assumed  available  single-target  MAP  tracker  to  compute  updated  parameters 

such  that 


jc(';'))e('"')|z 


:")} 


=  max 

. ^Tm 


(66) 


The  updated  state  parameters  are  MAP  estimates,  and  they  depend  on  the  current 

estimates  (64)  of  the  assignment  interference.  The  PMHT  algorithm  is  a  multitarget  tracking 
algorithm  that  uses  the  single-target  MAP  tracker  as  an  essential  building  block;  thus,  the  PMHT 
algorithm  is  conceptually  equivalent  to  an  iteratively  reweighted  bank  of  MAP  trackers. 
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The  computational  complexity  and  numerical  characteristics  of  the  PMHT  algorithm  are 
closely  related  to  the  complexity  and  numerical  accuracy  of  the  single-target  MAP  tracker.  The 
total  computational  effort  is  proportional  to  the  product  of  number  of  targets  M,  the  number  of 
PMHT  iterations  needed  to  estimate  the  assignment  interference  to  sufficient  accuracy,  and  the 
effort  needed  to  solve  the  single-target  MAP  parameter  estimates  from  (66).  Precise  estimates  of 
computational  effort  depend  on  the  particular  application. 

In  the  special  case  of  independent,  additive  white  Gaussian  state  and  measurement  noises, 
the  maximization  problem  (66)  is  equivalent  to  a  nonlinear  least-squares  problem.  This  is  easily 
seen  using  (45)  and  the  fact  that 

+  (67) 

is  equivalent  to 


^  (x  |x  ,  )  =  ■%  X  f  ,  {x  ,  ),Q  , 


(68) 


where  Q,  ,  is  the  covariance  of  the  white  Gaussian  noise  w,  ,  .  The  state  evolution  function 

^  and  the  a  priori  distribution  at  the  initial  batch  time  /  =  0  are  assumed  given.  Using  the 

Gaussian  a  priori  densities  (4),  the  centroids  (48),  and  the  effective  covariances  (49)  gives  the 
following  nonlinear  least-squares  problem  equivalent  to  (66): 


nun 

1*0, -^1. . *7>.} 


T  * 

(x  )l' (/?(';')  (x  )1 

tm  tnt\  tm/j  \  t-\,m  J  /  J 


(69) 


The  solution  of  (69)  is  the  updated  parameter  vector ^x|^'^'^,x|'J‘\...,x[,^'^j .  In  this  case,  a  special 

purpose  algorithm  such  as  Levenberg-Marquardt  (reference  10,  chapter  14.4)  can  be  used  to 
solve  (69)  and  can  serve  as  the  single-target  MAP  tracker.  Further  consideration  of  these  issues 
lies  outside  the  scope  of  this  report. 
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6.  THEORETICAL  ASPECTS  OF  THE  PMHT  ALGORITHM 


6.1  CONVERGENCE  PROPERTIES  OF  THE  PMHT  ALGORITHM 

Algorithms  derived  from  the  EM  method  have  two  useful  properties.  These  algorithms  are 
guaranteed  to  converge  under  very  mild  assumptions,  and  they  typically  get  very  close  to  the  limit 
point  during  their  first  few  iterations,  even  when  poorly  initialized.  It  is  known  that  the 
convergence  rate  of  EM  algorithms  is  only  asymptotically  linear;  however,  these  results  do  not 
explain  the  practical  observation  of  rapid  early-convergence  rate.  Considering  the  various 
modeling  approximations  that  are  made  in  practical  problems,  it  is  unlikely  that  exact  solutions  are 
really  necessary  for  the  multitarget  tracking  application  —  a  good  approximation  is  typically 
sufficient. 

If  greater  precision  is  required  than  is  provided  by  computing  a  few  steps  of  the  EM  method, 
it  may  be  wiser  to  use  a  compound  computational  procedure:  Begin  with  EM  to  get  close  to  the 
solution  quickly,  and  then  switch  to  a  more  rapidly  convergent  algorithm  (e  g.,  Newton-Raphson). 
The  choice  of  optimization  algorithm  is  theoretically  irrelevant,  as  long  as  ML  estimates  are 
computed  from  the  marginal  PDF  (22),  the  basis  of  the  PMHT  formulation.  The  EM  method  was 
chosen  in  this  study  for  this  estimation  problem  because  it  is  structurally  compatible  with  the 
notions  of  the  single-target  Kalman  filter. 

At  every  iteration,  the  PMHT  algorithm  is  guaranteed  to  either  increase  the  likelihood  of  the 
marginal  PDF  P(Z,X),  defined  by  equation  (22),  or  to  be  converged  to  a  stationary  point  of  the 
marginal  PDF.  This  property  alone  does  not  guarantee  convergence,  however.  For  example,  if 
the  marginal  is  unbounded  above,  the  PMHT  iterates  may  be  attracted  to  infinite  singularities  of 
the  likelihood  function,  and  this  can  happen  even  though  each  iterate  is  well  defined 
mathematically.  The  general  theory  (see  reference  6)  shows  that  boundedness  from  above 
guarantees  convergence  to  either  an  ML  point  or  a  stationary  point  of  the  marginal  PDF  (22). 
Although  stationary  points  may  seem  of  little  practical  importance,  convergence  to  local  ML 
solutions  instead  of  a  global  ML  solution  may  sometimes  cause  difficulties.  Section  9  discusses  a 
meaningful  local  solution  for  an  example  involving  two  crossing  targets.  In  a  multitarget  tracking 
application,  PMHT  requires  sliding  the  batch  as  new  measurement  scans  arrive,  consequently,  it  is 
probable  that  careful  initialization  of  the  starting  point  of  the  PMHT  algorithm  for  each  new 
batch,  using  the  solution  from  the  previous  batch,  will  result  in  convergence  to  global  ML 
solutions.  This  aspect  of  PMHT  requires  further  research. 

Under  linear-Gaussian  statistics,  an  appropriate  observability  condition  guarantees  that  the 
PMHT  marginal  PDF  is  bounded.  One  such  condition  is  that  the  individual  target-error- 
covariances  are  positive  definite  throughout  the  batch  for  all  possible  measurement-to-track 
assignments  K.  There  is,  however,  an  easily  overlooked  subtlety.  Consider  target  m  in  isolation. 

If  even  one  of  the  process  gain  matrices  |  (first  defined  in  equation  (1))  is  less  than  full  rank, 
then  the  PDF  (38)  is  unbounded  because  at  least  one  of  the  process  PDFs  „j|  is  a 

degenerate  Gaussian  distribution.  It  is  a  simple  matter  to  require  all  the  matrices  to  have 
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full  rank  to  avoid  the  degeneracy;  however,  most  kinematic  process  models  have  deficient  rank, 
so  this  requirement  is  overly  restrictive.  (Increases  in  the  marginal  PDF  are  easily  verified 
numerically  under  full-rank  conditions,  and  this  serves  in  practice  to  check  the  implementation  of 
the  PMHT  algorithm.)  Fortunately,  the  marginal  PDF  (22)  is  bounded  above  on  the  linear 
subspace  in  which  the  individual  target  processes  are  confined.  Because  the  PMHT  estimates  are 
themselves  confined  to  the  same  linear  subspace,  the  PMHT  algorithm  is  guaranteed  to  converge. 
Further  study  of  this  topic  is  outside  the  scope  of  this  report. 


6.2  FISHER  INFORMATION  MATRICES  FOR  PMHT 

The  PMHT  error-covariance  matrices,  or  inverse  FIMs,  for  the  target  states  appear 
incidentally  as  a  by-product  of  the  M-step  of  the  PMHT  algorithm,  and  it  is  unclear  from  the  EM 
method  context  exactly  how  these  matrices  should  be  interpreted.  The  difficulty  stems  from  the 
fact  that  the  cross-entropy  function  is  defined  for  its  desirable  analytical  properties,  and  not 
derived  from  statistical  assumptions  concerning  the  measurements.  Consequently,  quantities,  such 
as  FIMs,  that  are  computed  from  the  cross-entropy  function  lack  statistical  interpretation  in  the 
general  theory  of  the  EM  method.  The  purpose  of  this  subsection  is  to  present  one  such  statistical 
interpretation  for  the  FIMs  defined  for  the  PMHT  algorithm.  The  interpretation  presented  is  not 
theoretically  completely  satisfactory,  and,  therefore,  further  work  on  this  topic  is  justified. 

The  PDF  (40)  is  easily  interpreted  statistically  in  terms  of  a  Markov-state  sequence  and  the 
scan  PDF  .  The  difficulty  lies  in  the  statistical  interpretation  of 

is  not  equivalent  to  a  PDF  that  corresponds  to  making  hard  assignments  because  the  weights  (24) 
do  not  converge  to  hard  assignment  decisions;  that  is,  no  weight  converges  to  either  zero  or  one. 

The  interpretation  adopted  is  that  E  (z  lx  )  is  the  PDF  of  a  randomized  decision  rule  defined 

over  the  ensemble  of  all  possible  batch  measurements  {Z).  Each  member  of  the  ensemble  is 
assumed  to  have  exactly  the  same  number  of  measurements  in  every  scan  as  the  given  batch 
measurement  Z,  the  only  measurement  available  in  practice.  For  each  member  of  the  ensemble, 
hard  measurement-to-track  assignments  are  made,  however,  assignments  are  discrete  random 
variables,  and  the  probability  of  different  assignments  are  given  by  the  weights  (24).  The  PMHT 
algorithm  and  the  given  batch  measurement  Z  are  used  to  obtain  the  weights  (24),  and  these 
weights  are  applied  to  the  ensemble  (Z).  The  appropriate  PDF  for  a  randomized  decision  rule 

conceptualized  in  this  manner  is  the  scan  PDF  E^^Zjx^  j .  The  FIM  corresponding  to  this 
interpretation  of  the  measurements  is  computed  by  the  PMHT  algorithm  from  the  function  (40). 

FIMs  for  unbiased  estimators  are  derived  directly  from  a  given  joint  PDF.  The  joint  PDF 
may  incorporate  nonrandom  parameters,  random  parameters  (see  reference  1 1,  pp.  72,  84-5),  and 
mixed  random  and  nonrandom  parameters.  Unfortunately,  the  FIM  is  undefined  for  the  PMHT 
observer  PDF  (17)  because  the  full  parameter  list  comprises  the  random  parameters  X  and  K,  and 
the  gradient  with  respect  to  K  does  not  exist  because  K  is  discrete.  This  difficulty  is  circumvented 
by  deriving  the  FIM  from  the  observer  marginal  PDF  P{Z,X.Y\),  defined  by  (22),  by  treating  it  as 
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a  joint  PDF  whose  full  parameter  list  comprises  the  random  parameter  A"  and  the  nonrandom 
parameter  IT.  This  FIM,  called  the  marginal  FIM,  for  X  and  IT  is  defined  by 


J(n)  jj  Kn{v;„  \0iP(z,X:n)]y(z,x.-n)dZdx 

Z  X 


n=n 


(70) 


where  FI  is  the  true  value  of  IT .  The  marginal  FIM  for  X  and  FI  is  very  difficult  to  evaluate 
explicitly  in  the  multitarget  case  because  the  targets  are  coupled  by  assignment  interference,  as 
discussed  in  section  2.  A  study  of  the  marginal  FIM  is  outside  the  scope  of  this  report. 
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7.  VARIATIONS  AND  EXTENSIONS  OF  PMHT 


7.1  TARGET  MEASUREMENT  PROBABILITIES,  n 

Several  interesting  variations  of  the  PMHT  likelihood  function  have  cross-entropy  0 
functions  that  are  easily  obtained.  Because  of  the  decoupling  of  the  estimation  steps  for  IT  and  X, 
these  variations  result  in  only  slight  changes  in  the  PMHT  algorithm  and  are  readily  implemented. 
This  subsection  discusses  two  alternatives  for  11. 


Some  applications  may  require  that  the  target  measurement  probability  vectors  be 

identical  across  the  batch,  but  not  assumed  known  a  priori.  The  BIN  for  this  variation  is 
unchanged  from  that  of  figure  1  because  it  is  equivalent  to  a  constraint  on  the  parameterization  of 
the  joint  PDF;  that  is,  the  conditional  independence  assumptions  of  the  random  variables  are 

unaltered  by  this  constraint,  but  are  only  reparameterized.  In  this  case,  the  [0,  n}  terms  (29)  are 
modified  by  substituting 


n  =n,  -n.  = 

j  2j 


= 


s  = 


(71) 


The  probability  n ^  is  used  only  in  this  section,  and  it  is  not  to  be  confused  with  the  probability 
vector  TT^  given  by  equation  (9).  The  constraint  (71)  modifies  the  M-step  slightly.  Using  the 

notational  conventions  of  section  5,  the  updating  recursions  (46)  through  (49)  are  replaced  by  the 
recursions 


w 


s=\ 


(72) 
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R 


(75) 


respectively,  where  the  constant  V  in  (73)  is  given  by 
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(76) 


N  =  n^+  ■■■  +  72^ . 

The  recursions  (50)  through  (57)  are  unchanged.  The  appropriate  equations  in  the  general  case 
are  derived  similarly.  This  variation  is  not  considered  further. 

If  the  target  measurement  probabilities  11  are  assumed  known  a  priori,  for  example, 

7c^^  =  \l  M  for  all  /  and  m,  then  the  Q  function  for  this  case  comprises  the  terms  |  iri 

equation  (30),  but  it  does  not  include  the  terms  n }  M-step  in  this  case  updates  X  in 

exactly  the  same  way  as  before,  but  it  uses  the  specified  a  priori  n  instead  of  the  estimated  ft- 
More  generally,  the  0  function  is  easily  modified  to  incorporate  a  specified  a  priori  distribution 
on  n,  if  one  is  available.  These  variations  are  not  considered  further  in  this  report. 


7.2  ADAPTIVE  COVARIANCE  ESTIMATION 


The  measurement  and  target-covariance  matrices  and  are  assumed  known  in 

the  above  development  for  the  linear-Gaussian  case.  When  these  matrices  are  unknown,  they  can 
be  treated  as  parameters  to  be  estimated  from  the  batch  measurement  data.  Although  covariance 
estimation  problems  are  outside  the  scope  of  this  report,  the  necessary  conditions  are  easily 
obtained  and  are  presented  below  because  of  their  potential  use. 

The  EM  approach  to  covariance  estimation  does  not  alter  the  likelihood  structure  (15)  or  the 
cross-entropy  Q  function,  but  it  does  significantly  alter  the  M-step  because  the  covariance  matrix 
estimates  are  coupled  to  the  state  estimation  equations.  The  0  function  for  this  case  is  identical 
to  (28),  but  its  gradient  with  respect  to  each  of  the  matrices  R  and  0  must  now  be  taken  and 

tn7 

set  to  0.  Taking  the  gradient  of  (30)  with  respect  to  applying  the  identity  (reference  5, 
section  2.14.2) 


V,  \o^7t{x\p,T)  =  -S-'  +  Z-'(x-/i)(x-/2)'Z-', 


and  setting  the  result  to  zero  gives  the  estimate 


R  = 

tm 


f 

s 

Vr=l 


(77) 


The  estimate  (77)  is  full  rank  (with  probability  one)  if  the  dimension  N ^  of  the  measurements  is 
less  than  the  number  of  measurements.  Similarly,  assuming  the  matrices  are  nonsingular,  the 
gradient  of  (30)  with  respect  to  gives  the  estimate 
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(78) 


The  estimate  (78)  is  fiill  rank  only  if  the  dimension  of  the  target-state  variable  N  is  exactly  1. 

Rank  difficulties  for  covariance  matrix  estimation  within  a  scan  are  significantly  reduced  by 
requiring  the  covariance  matrices  to  be  stationary  throughout  the  batch,  that  is,  and 

Qtm  =  Qm  ^  additional  requirement  gives  the  measurement  covariance  estimate 
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(79) 


The  covariance  matrix  estimate  (79)  is  theoretically  of  full  rank  (with  probability  one)  only  if  the 
the  total  number  iV  of  individual  measurements  is  greater  than  the  dimension  of  the 

measurement  variables.  Similarly,  assuming  =  G^,  the  estimate  for  the  process  covariance  is 
given  by 

1  r  ^ 

Q  =lyfG‘’(f  ,  -F  X  )1[G''(x  ,  -F  X  )1  .  (80) 

The  estimate  (80)  is  the  average  of  the  estimates  (78),  and  it  is  theoretically  of  full  rank  (with 
probability  one)  if  Fis  greater  than  the  dimension  N of  the  state  variables  for  target  m. 

Equations  (77)  through  (80)  are  nonlinearly  coupled  with  the  M-step  state  estimates,  but 
they  are  decoupled  into  M  smaller  systems,  one  for  each  target  model.  An  explicit  solution  of 
these  nonlinear  equations  is  not  available.  However,  it  is  straightforward  to  show  that  the 
estimates  (79)  and  (80)  are  the  basic  equations  for  the  generalized  EM  algorithm  (reference  6)  for 
estimating  covariances. 

Implementing  the  covariance  estimates  in  the  mathematical  forms  given  above  will 
unnecessarily  square  the  numerical  condition  number  of  the  underlying  "data  matrix."  Numerical 
ill-conditioning  arising  from  this  source  is  completely  avoided  by  the  use  of  QR  methods 
(reference  12).  QR  algorithms  are  stable  numerically  and  require  very  little  additional 
computational  effort,  so  they  are  highly  recommended  for  avoiding  sample  covariance  matrix 
formation.  Another  source  of  numerical  ill-conditioning  for  the  estimate  (79)  is  the  reduction  of 
the  effective  numerical  rank  that  occurs  when  one  or  more  of  the  summand  coefficients  becomes 
extremely  small.  QR  algorithms  can  limit,  but  not  fully  overcome,  ill-conditioning  due  to  the 
dynamic  range  of  these  coefficients. 
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8.  COMPARISONS  WITH  OTHER  METHODS 


8.1  MAXIMUM  LIKELIHOOD  DATA  ASSOCIATION 

Avitzour,  in  reference  13,  considers  data  association  using  an  ML  algorithm  that  is  based  on 
the  EM  method.  Several  of  the  ideas  described  are  closely  related  to  the  ideas  developed  in  this 
report;  however,  there  are  also  significant  differences.  The  points  of  similarity  are  discussed  first. 
Both  methods  consider  data  association  for  multitarget  tracking  fi-om  a  probabilistic  perspective, 
and  both  consider  batch  measurements.  Both  also  use  the  EM  method  to  estimate  target  states 
and  assignment  probabilities.  Moreover,  both  recognize  the  significance  in  the  multitarget 
tracking  application  of  the  assignment  probabilities  (24).  Avitzour  is  “not  concerned  with  finding 
the  correct  association  of  measurements  to  targets/clutter”  because  assignment  probabilities  are 
“obtained  as  a  by-product”  of  state  estimation.  This  view  of  measurement  assignments  is  identical 
to  the  view  expressed  in  this  report.  Finally,  both  discuss  data  association  using  hard  assignments 
as  joint  MAP  estimation  of  states  and  discrete  random  variables  called  assignments  and  derive 
MAP  state  estimates  from  a  marginal  distribution  over  the  assignments  using  the  EM  method. 

One  of  the  important  differences  between  reference  13  and  this  report  is  that  Avitzour  does 
not  fully  interpret  the  EM  method  in  the  context  of  the  multitarget  tracking  application.  In 
particular,  an  observer  whose  state  includes  the  assignments  is  not  defined  explicitly  in  reference 
13.  It  is  this  important  feature  of  the  observer  that  distinguishes  PMHT  from  joint  probabilistic 
data  association  (JPDA)  and  MHT,  enables  the  EM  method  to  be  applied  effectively,  and  provides 
a  unified  cogent  exposition  of  the  probabilistic  issues. 

Error-covariance  estimates  are  estimated  and  discussed  in  this  report,  but  are  absent  from 
reference  1 3 .  The  primary  reason  this  study  can  estimate  covariances  is  the  recognition  that  the 
exponential  of  the  cross-entropy  function  is  proportional  to  the  likelihood  structure  of  a  Kalman 
filter  (compare  with  equations  (38)  and  (40))  that  treats  an  entire  measurement  scan  2t  as  a  single 
measurement.  This  important  relationship  between  the  M-step  and  the  Kalman  filter  is  not 
recognized  in  reference  1 3 . 

The  central  role  played  by  conditional  independence  in  multivariate  statistical  problems  such 
as  multitarget  tracking  is  not  emphasized  by  Avitzour.  The  BIN  graphical  representation  of 
conditional  independence  structure  of  the  PMHT  approach  is  novel  to  this  study  and  is  not  used  in 
reference  13.  The  BIN  graph  helps  clarify  the  differences  between  PMHT  and  JPDA. 

In  reference  13,  Avitzour  does  not  allow  process  noise,  although  he  does  allow  moving 
targets.  This  is  not  a  serious  difference,  however,  as  it  is  possible  to  extend  his  method  to  include 
process  noise.  Finally,  Avitzour  presents  a  potentially  useful  termination  criterion  for  the  EM 
method.  This  termination  criterion  may  be  useful  in  this  study  also,  but  it  is  not  used  here. 
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8.2  JOINT  PROBABILISTIC  DATA  ASSOCIATION 


In  the  case  of  linear-Gaussian  statistics,  the  PMHT  algorithm  resembles  the  JPDA 
algorithm;  however,  PMHT  is  fundamentally  different  from  JPDA.  The  resemblance  is  closest 
for  a  unit  batch  of  size  {T=  1),  which  is  the  batch  size  assumed  for  PMHT  in  the  present 
discussion.  Linear-Gaussian  statistics  are  assumed  thoughout  this  discussion. 

Both  PMHT  and  JPDA  assume  that  the  measurements  within  a  scan  are  independent 
samples  of  a  random  variable  with  a  mixture  Gaussian  PDF.  JPDA  assumes  a  priori  that  the 
scan  measurement  mixture  is  uniform-Gaussian,  that  is,  the  mixing  proportions  of  the  Gaussian 
components  are  equal  and  their  means  are  the  predicted  measurement  means.  JPDA  use  of 
predicted  states  is  one  source  of  track  bias.  In  contrast  to  JPDA,  PMHT  uses  ML  estimates  of 
the  mixing  proportions  and  the  means  of  the  Gaussian  components  in  the  scan  measurement 
mixture.  Because  targets  may  not  be  equally  represented  by  the  measurements,  PMHT  use  of 
nonuniform  mixing  proportions  permits  improved  modeling  of  scan  measurements.  The  ML 
estimates  of  the  component  means  are  the  estimated  target  measurement  means  for  the  current 
scan.  Because  PMHT  uses  estimated  -  as  opposed  to  predicted  (as  in  JPDA)  -  target 
measurement  means,  PMHT  target  state  estimates  should  often  have  smaller  biases,  smaller 
critical  separation  distances  at  which  tracks  coalesce,  and  better  tolerance  to  target  maneuvers 
than  JPDA. 

JPDA  first  tracks  measurements  and  then  smoothes  the  resulting  state  estimates  using 
convex  combinations.  The  convex  combinations  yield  the  output  state  estimates  and  their 
corresponding  error  covariances.  Convex  combinations  occur  in  the  JPDA  state  space  because 
of  the  way  in  which  JPDA  treats  assignments  as  events.  In  contrast,  PMHT  first  smoothes 
measurements  using  convex  combinations,  and  then  tracks  the  smoothed  measurement.  The 
convex  combinations  of  measurements  are  the  probabilistic  measurement  centroids  (62),  and  the 
centroids  constitute  synthetic  measurements  to  be  tracked.  Convex  combinations  occur  in  the 
PMHT  measurement  space  because  assignments  are  incorporated  into  the  observer  (compare 
with  equation  (22)). 

PMHT  error  covariances  are  different  from  those  of  JPDA.  Covariances  computed  by 
JPDA  are  readily  interpreted  statistically  because  JPDA  treats  assignments  as  events. 
Unfortunately,  a  similar  interpretation  is  not  applicable  to  PMHT  covariances.  The  PMHT 
covariances  are  tied  implicitly  to  the  use  of  the  EM  method  to  compute  ML  state  estimates  from 
the  marginal  PDF  (22).  For  further  discussion  see  section  6.2. 

Measurement  gates  are  used  to  censor  measurements  that  are  too  far  from  predicted 
measurements  to  be  associated  with  targets.  Gates  will  probably  always  be  needed  to  protect 
theoretical  mathematical  models  from  the  vagaries  of  real  data  as  well  as  to  improve  computation 
time  in  practical  systems.  JPDA  is  explicitly  formulated  using  gates,  but  PMHT  does  not  use 
gating  in  its  formulation.  However,  gates  are  readily  derived  for  PMHT  by  thresholding  the 
assignment  weights  (24). 
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Gates  are  equivalent  to  applying  zero  weights  to  certain  measurements.  JPDA  and  PMHT 
discount  zero-weighted  measurements  in  very  different  ways.  In  JPDA,  the  convex  combination 
of  state  estimates  is  unaffected  by  the  addition  of  zero-weighted  estimates.  In  PMHT,  zero- 
weighted  measurements  contribute  a  multiplicative  factor  of +1  to  the  conditional  PDF  defined  by 
equation  (41).  Consequently,  zero-weighted  measurements  do  not  affect  PMHT  state  estimates. 
The  robustness  of  the  PMHT  and  JPDA  state  estimates  to  measurements  with  small  weights  (i.e., 
outliers)  remains  to  be  studied. 


8.3  GAUSSIAN  SUM  TRACKING  FILTERS 

Sorenson  and  Alspach  (references  14  and  15)  describe  a  Gaussian  sum  tracking  filter  for 
nonlinear  estimation.  This  approach  is  also  discussed  in  reference  5,  section  8.4.  The  problem 
they  studied  was  restricted  to  single  targets  -  multiple  targets  are  not  mentioned.  The  main 
reason  to  mention  multiple  targets  in  this  report  is  that  the  PMHT  marginal  PDF  defined  by 
equation  (22)  can  be  thought  of  as  a  Gaussian  sum  approximation  of  a  nonlinear  measurement 
PDF.  This  interpretation  suggests  that  the  PMHT  approach  may  be  applied  to  general  nonlinear 
estimation  problems. 

The  general  nonlinear  estimation  problem  is  investigated  in  references  14  and  15  as  a 
problem  in  the  approximation  of  the  posterior  PDF.  The  form  of  PDF  approximation  used  is  that 
of  a  Gaussian  mixture,  so  this  approach  is  fundamentally  different  from  the  extended  Kalman 
filter.  In  reference  14,  the  main  problem  discussed  is  that  of  a  system  with  linear  process  and 
measurement  equations  with  additive  non-Gaussian  noises.  A  Gaussian  mixture  approximates  the 
prior  density,  and  the  additive  process  and  measurement  noises  have  PDFs  that  are  approximated 
by  Gaussian  mixtures.  Under  these  assumptions,  the  posterior  density  is  shown  to  be  exactly  a 
Gaussian  mixture.  In  reference  15,  nonlinear  process  models  are  approximated  by  Gaussian 
mixtures  by  linearizing  the  process  model  about  the  mean  of  each  Gaussian  in  the  mixture.  The 
primary  conceptual  drawback  to  this  approach  is  that  the  number  of  Gaussian  components  in  the 
mixture  representing  the  posterior  PDF  grows  rapidly  as  time  evolves;  therefore,  pruning  terms 
from  the  Gaussian  mixture  is  necessary  in  practice. 

Sorenson  and  Alspach  estimated  the  parameters  of  the  approximating  mixtures  by  applying 
gradient-based  numerical  methods  because,  at  the  time  their  papers  were  written,  the  method  of 
EM  had  not  been  developed  in  generality  and  was  largely  unknown  outside  the  statistical 
community.  The  primary  similarity  with  PMHT  lies  in  their  recognition  that  the  state  estimate 
produced  by  the  Gaussian  sum  tracking  filter  is  obtained  as  a  convex  combination  of  several 
linear-Gaussian  filters.  The  filters  Sorenson  and  Alspach  refer  to,  however,  are  independent  filters 
that  are  not  coupled  by  assignment  interference. 


8.4  CENTROID  GROUP  TRACKING 

Blackman,  in  section  1 1.2  of  reference  16,  discusses  tracking  groups  of  closely  spaced 
targets  by  tracking  the  centroid  of  the  group.  This  application  is  a  specialized  multitarget  tracking 
problem,  and  the  approach  described  is  quite  different  from  PMHT.  Nonetheless,  there  are 
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similarities  in  certain  details.  In  particular,  the  centroids  of  the  measurements  are  used  as 
synthetic  measurements  for  updating  the  group  track  centroid.  It  is  necessary  in  this  approach  to 
gate  the  measurements  before  forming  the  measurement  centroids,  and  such  gating  is  tantamount 
to  making  hard  assignments  of  measurements  to  different  groups.  The  PMHT  equivalent  of  hard 
assignments  and  gates  is  to  threshold  the  PMHT  weights  (24).  There  are  also  similarities  between 
determining  the  number  of  targets  in  a  group  and  the  choice  of  the  parameter  M  in  the  PMHT 
algorithm,  a  model-order  sdection  problem  that  falls  outside  the  scope  of  this  report. 

Formation  group  tracking  is  also  discussed  by  Blackman  in  reference  16,  section  1 1.3.  In 
formation  tracking  the  targets  are  not  assumed  to  be  independent,  and  the  availability  of 
information  concerning  their  coordinated  behavior  is  exploited  to  advantage.  Although  the 
PMHT  approach  is  not  limited  to  independent  targets,  this  topic  also  falls  outside  the  scope  of  this 
report. 
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9.  CROSSING  TARGETS  WITH  LINEAR-GAUSSIAN  STATISTICS 


Section  9  provides  an  example  that  demonstrates  PMHT  ability  to  estimate  measurement- 
to-track  assignments.  It  is  given  that  two  targets  (M=  2)  are  present  and  constrained  to  move  in 
the  xy  plane.  For  simplicity,  positions  are  given  in  meters.  Both  target  tracks  begin  at  /  =  0,  are 
simulated  without  process  noise,  and  have  constant  velocity.  Target  1  is  heading  45°  from  the 
+x-axis,  has  speed  +1  m/s,  and  begins  its  track  at  the  origin  (0,  0)  of  the  x-y-plane.  Target  2  is 

heading  +90°  from  the  +x-axis,  has  speed  V2  /  2  »  0.707  m/s,  and  its  track  begins  at  +12.5  m  on 
thex-axis.  Given  this  geometry,  the  targets  become  superposed  at  12.5-\/2  «  17.68  s. 
Measurement  scans  are  taken  at  1 -second  intervals,  and  the  batch  length  is  T=  25  s.  One 
measurement  is  simulated  from  each  target  for  t  =  1,2,... ,25  s.  No  measurements  are  simulated  at 
/  =  0.  For  simplicity,  false  alarms  are  not  simulated,  and  each  target  generates  exactly  one 
measurement  in  the  simulation.  The  number  of  measurements  per  scan  is  always  =  2,  so  that 

the  total  number  of  measurements  in  the  batch  is  50.  The  state  vector  of  each  target  comprises 
the  x-  and  y-components  of  position  and  velocity.  For  estimation,  the  PMHT  target  state  models 
are 


m,f+l 

"1 

1 

0 

o’ 

mt 

X 

m./  +  l 

_ 

0 

1 

0 

0 

-f  cr^ 

state 

mt 

•^m.r+l 

0 

0 

1 

1 

mt 

0 

0 

0 

1 

1 

1 _ 

- 1 

m=l,2. 


(81) 


where  the  process  noise  is  white  Gaussian  with  covariance  matrix  given  by  the  4  x  4-identity 
matrix  scaled  by  the  dimensionless  quantity  =  10“^ .  The  Gaussian  a  priori  distributions 
(compare  with  equation  (4))  for  targets  1  and  2  have  mean  vectors 

Target  1 :  (0  m,  a  m/s,  0  m,  a  m/s) 

Target  2;  (12.5  m,  0  m/s,  0  m,  a  m/s), 

where  the  dimensionless  constant  a  =  V2  /  2  .  The  covariances  of  the  a  priori  distributions  are 
equal  to  the  4  x  4-diagonal  matrix  c//ag|l6  m\  9  m*  /  s\  16  m\  9  mV  s^  j  for  both  targets. 

Prior  distributions  are  centered  on  the  simulated  target  state  at  /  =  0,  but  they  have  large 
uncertainties  in  position  and  velocity.  Because  the  simulation  is  defined  in  such  a  way  that  each 
target  generates  one  and  only  one  measurement,  the  measurement  equations  can  be  written  in  the 
form 
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(82) 
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where  the  measurement  noise  is  white  Gaussian  with  covariance  matrix  given  by  the  2  x  2-identity 
matrix  scaled  by  the  dimensionless  quantity  .  Equation  (82)  is  obtained  by  substituting  r  =  m 
in  equation  (7);  however,  the  target  labeling  implicit  in  (82)  is  not  used  by  the  PMHT  algorithm. 

It  is  unnecessary  to  randomly  shuffle  the  measurements  {^,,,22)  before  calling  the 

PMHT  algorithm  because  data  order  does  not  alter  the  algorithm,  as  can  be  seen  from  equations 
(46)  through  (48). 

The  measurement  scale  cr]  in  equation  (82)  is  chosen  to  control  the  “effective”  duration  of 

the  track-crossing  time.  The  duration  of  the  track-crossing  event  is  defined  to  be  the  total  time 
the  two  targets  lie  within  three  measurement  noise  standard  deviations  in  the  x-y-plane.  Given  the 
above  geometry,  it  is  readily  seen  that  the  crossing  duration,  denoted  by  4,  is  linearly  proportional 
to  the  measurement  standard  error;  specifically,  the  cr^  » 8.485  The  measurement 

scalecr^  is  selected  so  that  tc  =  2  s,  hence,  cr^  =  0.2357.' 


Targets  1  and  2  are  initialized  in  the  correct  states  for  all  time  /  =  0, 1,...,25.  This  no-target- 
error  initialization  is  a  best-case  scenario  for  determining  the  number  of  iterations  to  convergence. 
Using  an  absolute  error  convergence  criterion  on  the  sequence  of  overall  likelihood  function 
iterates,  the  PMHT  algorithm  converged  in  38  iterations  when  the  convergence  error  s  =  0.0001. 
The  final  value  of  the  log  likelihood,  loge  P(Z,X),  was  +392.92.  Experience  with  the  PMHT 
algorithm  reveals  that  it  is  very  robust  to  poor  target  initialization;  therefore,  perfect  initialization 
is  unimportant  in  this  example.  (With  clutter,  however,  sensitivity  to  poor  target  initialization  is 
an  important  topic.  The  form  of  the  PMHT  algorithm  presented  here  assumes  no  clutter  in  order 
to  focus  clearly  and  solely  on  the  assignment  problem.  Incorporating  clutter  models  into  the 
PMHT  approach  is  straightforward  and  is  the  subject  of  an  ongoing  investigation.) 

Figure  3a  depicts  the  true  track  positions  overlaid  with  the  measurements,  and  figure  3b 
depicts  the  true  track  positions  overlaid  on  the  PMHT -estimated  track  positions.  Figure  4  depicts 
the  root-mean-square  error  in  the  estimated  track  position  and  speed  components  for  each  target. 
The  PMHT  assignment  weights  are  given  in  figure  5  for  both  targets  and  each  measurement.  The 
estimated  weights  are  clearly  indicative  of  the  correct  assignments  (known  only  via  simulation) 
everywhere  except  in  the  crossover  regime.  During  crossover,  the  evidence  from  the  data  is 
inconclusive  concerning  correct  assignments,  and  this  uncertainty  is  reflected  by  the  nearly  equal 
weights  given  measurements  in  this  regime. 
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METERS/SECOND 


0  S  10  15  20  0  5  10  15  20 

METERS  METERS 

a.  Simulated  Target  Tracks  and  Measurements  b.  Simulated  and  Estimated  Target  Tracks 


Figure  3.  Two  Target  Crossing  Track  Example  with  U  =  2 


SECONDS  SECONDS 

a.  Target  1  Position  Error  (T op)  b.  Target  2  Position  Error  (Top) 

and  Speed  Error  (Bottom)  and  Speed  Error  (Bottom) 

Figure  4,  PMHT  Position  and  Speed  Estimation  Errors  for  Each  Target 
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b.  Weights  for  Target  2  with  Zti  (Top) 
and  Zt2  (Bottom) 


Figure  5.  Converged  PMHT  Assignment  Weights  (18)  for  Two  Measurements  (82) 


Poor  PMHT  initialization  may  result  in  convergence  to  suboptimal  solutions;  i.e.,  the  PMHT 
algorithm  may  converge  to  only  a  local  ML  solution  -  not  the  global  solution  -  to  the  tracking 
problem.  In  the  multitarget  tracking  application,  these  suboptimal  solutions  may  be  meaningful 
because  more  than  one  interpretation  of  the  measurements  is  often  possible.  A  multiplicity  of 
meaningful  local  solutions  appears  to  be  an  important  feature  of  good  probabilistic  models  of  the 
problem;  If  the  observer  PDF  has  only  one  peak,  namely  the  global  peak,  then  it  does  not  contain 
within  it  alternative  interpretations  of  the  measurements.  Alternative  interpretations  are  conceived 
as  local  ML  solutions  whose  likelihoods,  determined  by  the  PMHT  algorithm,  are  significantly 
lower  than  the  likelihood  of  the  global  solution.  This  view  is  discussed  quantitatively  for  the  two 
target-crossing  examples  presented  above. 

Two  ML  solutions,  one  global  and  the  other  local,  are  given  in  table  1 .  These  solutions 
were  obtained  using  identical  a  priori  distributions  (for  all  values  of  L)  and  measurement  sets  (for 
each  value  of  4),  but  different  initializations  of  the  PMHT  algorithm.  The  initializations  used  to 
obtain  these  solutions  were  obtained  by  swapping  portions  of  the  true  tracks  (known  from 
simulation)  in  the  obvious  manner.  The  solution  likelihoods  given  in  table  1  are  estimated  using 
the  PMHT  algorithm  and  are  values  of  the  marginal  PDF  (22).  The  number  of  iterations  required 
for  PMHT  convergence  is  also  given  in  table  1 .  The  three  columns  in  table  1  correspond  to  easy 
(4  =  1),  moderately  difficult  (4  =  4),  and  difficult  (4  =  8)  crossing-trajectory  problems,  where 
problem  difficulty  is  quantified  by  the  crossing-time  duration  4.  As  seen  from  table  1,  for  a  given 
level  of  difficulty,  the  solutions,  ranked  by  decreasing  likelihood,  are  in  the  same  rank  order  as 
was  anticipated  intuitively.  The  likelihood  ratio  shows  that  the  crossing-track  solution  is 
significantly  more  likely  than  the  switching-track  solution.  The  likelihood  ratio  also  shows  that 
the  dynamic  range  of  the  likelihoods  of  the  two  solutions  decreases  with  increasing  problem 
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difficulty.  Interpreting  solution  likelihood  as  solution  explanatory  power,  the  decrease  in  dynamic 
range  means  that  the  local  solutions  provide  less  differentiated  explanations  of  the  data  as  the 
tracking  problem  becomes  more  difficult.  In  the  limit  as  >  qo  ,  all  local  solutions  provide 

equally  good  explanations  of  the  data. 


Table  1.  Local  ML  Tracking  Solutions  Ranked  by  Total  Likelihood  and 

Crossing-Track  Duration,  tc 


Solution 

Loge  P{Z,X)  and  Number  of  PMHT  Iterations 

Track  description 

^||||||||||^^ 

L=8 

(difficult) 

Crossing  tracks 

469.20 

3  iterations 

327.54 

29  iterations 

258.22 

34  iterations 

Switching  tracks 

459.62 

3  iterations 

318.88 

52  iterations 

252.88 

30  iterations 

Likelihood  ratio 

1.45x10'' 

5.77x10' 

2.09x10' 

The  solutions  listed  in  table  1  are  not  the  only  solutions  to  the  problem.  In  particular, 
because  the  observer  PDF  is  unchanged  by  interchanging  the  roles  of  targets  1  and  2,  two  other 
solutions  are  readily  obtained.  These  alternative  solutions,  however,  cannot  be  considered 
different  from  those  listed  in  table  1 .  The  interchangeability  of  certain  observer  PDF  parameters  is 
related  to  the  more  general  issue  of  parameter  “identifiability,”  as  it  is  known  in  the  general 
statistical  literature.  This  statistical  terminology  is  unrelated  to  the  system  identifiability 
terminology  in  the  Kalman  filtering  literature. 
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10.  CONCLUSIONS  AND  RECOMMENDATIONS 


10.1  CONCLUSIONS 

The  study  presented  in  this  report  is  a  probabilistic  approach  to  the  measurement-to-track 
assignment  problem.  Measurements  were  not  assigned  to  tracks  as  in  traditional  multi-hypothesis 
tracking  (MHT)  algorithms;  instead,  the  probability  that  each  measurement  belongs  to  each  track 
was  estimated  using  a  maximum  a  posteriori  (MAP)  method.  These  measurement-to-track 
probability  estimates  are  intrinsic  to  the  multitarget  tracker  called  the  probabilistic  multi¬ 
hypothesis  tracking  (PMHT)  algorithm.  The  PMHT  algorithm  is  computationally  practical 
because  it  requires  neither  enumeration  of  measurement-to-track  assignments  nor  pruning.  The 
PMHT  algorithm  is  an  optimal  MAP  multitarget  tracking  algorithm. 


10.2  TOPICS  OF  PRACTICAL  IMPORTANCE 

The  following  topics  of  practical  importance  related  to  the  PMHT  approach  remain  to  be 
investigated. 

1 .  The  choice  of  the  number  of  target  models  M  is  an  important  model-order  selection 
problem  that  is  of  practical  importance  to  a  mature  multitarget  tracking  algorithm.  It  is 
anticipated  that  PMHT  will  be  robust  to  mismatch  between  the  true  number  of  targets  and  the 
number  M  assumed  by  PMHT,  provided  that  M  is  at  least  as  great  as  the  true  number  of  targets. 
The  reason  for  suspecting  such  tolerance  for  PMHT  is  that  "extra"  target  models  can  be  used  to 
model  background  clutter  and  noise  levels. 

2.  Outputs  from  the  previous  batch  can  be  used  to  initialize  the  PMHT  algorithm  for  the 
current  batch,  but  this  procedure  does  not  reduce  computational  complexity  because  it  is  not 
recursive.  Development  of  a  recursive  PMHT  algorithm  that  operates  between  successive 
measurement  batches  is  an  important  subject  for  future  work. 

3.  The  robustness  of  PMHT  algorithms  against  target  coalescence  is  of  particular 
importance.  All  tracking  algorithms  will  coalesce  sufficiently  closely  spaced  targets.  Because 
PMHT  estimates  are  optimal  empirical  Bayesian  estimates,  it  is  anticipated  that  PMHT  ability  to 
resolve  closely  spaced  targets  will  be  at  least  as  good  as  other  currently  available  methods.  (A 
closely  related  topic  is  track  estimation  bias.) 

4.  The  sensitivity  of  the  PMHT  algorithm  to  target  maneuvers  must  also  be  studied. 

5.  Bayesian  priors  for  the  assignments  K  can  be  used  in  the  PMHT  approach,  if  they  are 
available.  Such  priors  might  be  found  by  exploiting  additional  target  information  that  is  not 
resident  in  the  sensor  measurements  Z;  however,  such  priors  have  not  been  studied. 
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10.3  TOPICS  OF  THEORETICAL  INTEREST 


The  following  topics  of  theoretical  interest  should  be  studied  further. 

1.  The  sequence  of  error  covariances  generated  by  the  PMHT  algorithm  is  significant 
because  multitarget  FIMs  are  not  readily  definable  using  other  approaches.  (See  section  6.) 

2.  Convergence  of  the  PMHT  algorithm  is  related  to  target  observability.  The  multitarget 
observability  condition  for  PMHT  stated  in  section  6  essentially  requires  that  all  targets  be 
individually  observable.  Conditions  under  which  it  may  be  possible  to  weaken  this  requirement 
would  be  very  interesting.  For  instance,  if  the  targets  are  not  independent  because  they  move  in  a 
coordinated  formation,  it  may  be  possible  to  improve  the  effective  target  observability. 

3.  The  relationship  between  the  linear-Gaussian  PMHT  algorithm  and  Gaussian  sum 
approximation  methods  (see  references  14  and  15)  for  single-target  tracking  problems  with 
nonlinear  non-Gaussian  state  and  measurement  processes  deserves  investigation.  In  this  case, 
PMHT  approximates  the  exact  nonlinear  state  distributions  by  a  (convex)  linear  superposition  of 
several  linear-Gauss-Markov  processes  that  are  linked  by  assignment  interference.  The  PMHT 
approach  may  result  in  efficient  algorithms  for  these  methods. 


10.4  TOPICS  OF  BOTH  THEORETICAL  AND  PRACTICAL  INTEREST 

The  following  topics  of  both  theoretical  and  practical  interest  merit  further  study. 

1 .  The  adaptive  estimation  of  target  covariances  (i  .e.,  )  and  the  measurement 

covariances  (i.e.,  using  the  EM  method  need  further  study.  The  equations  for  adaptive  ML 

covariance  estimates  are  given  in  section  7.  These  adaptive  equations  are  coupled  with  the 
PMHT  target-state  estimates  because  of  assignment  interference.  Algorithms  for  solving  these 
coupled  equations  efficiently  remain  to  be  investigated. 

2.  Non-EM  methods  for  solving  for  the  parameters  of  the  marginal  density  (22)  may  be  very 
useful  in  practice,  but  have  not  been  considered  in  this  report. 
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